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ABSTRACT 

The  Method  of  Partial  Images  has  been  successfully 
applied  to  electrostatic  problems  involving  conductors  on 
a  dielectric  substrate.   This  same  method  is  investigated 
for  its  adaptation  to  dynamic  problems'.   A  Green's  Function 
is  derived  and  applied  to  the  problem  of  a  wire  dipole 
antenna  on  a  dielectric  substrate.   The  input  admittance 
of  the  antenna  is  computed  by  the  Method  of  Moments. 
Experimentally  measured  values  of  input  admittance  are 
compared  with  the  theoretical  values  and  the  error  is 
discussed. 
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I.   INTRODUCTION 

If  a  monopole  antenna  is  mounted  on  a  perfectly 
conducting  ground  plane,  it  will  reflect  totally  above  the 
plane  and  will  produce  an  image  below  the  plane.   This 
image  will  have  an  image  current  in  the  same  direction 
as  the  monopole  in  order  to  satisfy  boundary  conditions 
at  the  ground  plane.   In  essence,  because  of  its  image, 
the  monopole  will  act  as  a  dipole  in  a  free-space.   Hence, 
if  the  ground  plane  is  removed,  the  image  may  be  replaced 
by  another  physically  real  monopole  thereby  creating  a  dipole 
In  other  words,  conditions  above  the  reflective  plane  remain 
unchanged  if  the  plane  is  removed  and  a  real  object  replaces 
the  image. 

A.   THE  USE  OF  IMAGE  THEORY 

With  the  monopole,  the  thickness  of  the  ground  plane 
is  assumed  to  negligible.   This  results  in  only  one  image.  . 
However,  if  the  ground  plane  is  a  dielectric  with  finite 
thickness,  total  reflection  no  longer  occurs.   Instead, 
images  are  produced  from  partial  reflection  at  the  surface, 
reflection  within  the  dielectric  and  transmission  through 
the  dielectric.   These  images  are  used  in  the  "method  of 
partial  images"  [1].   When  the  dielectric  region  is  removed, 
these  images  may  be  considered  real  charges  and  are  summed 
as  an  infinite  series.   However,  a  finite  number  of  them 
give  good  results. 


B.   A  DYNAMIC  SOLUTION 

The  method  of  partial  images  has  been  successfully 
applied  to  static  problems.   The  purpose  of  this  research 
is  to  investigate  its  application  in  solving  dynamic 
problems.   A  search  for  a  Green's  Function  is  conducted. 
Once  found,  it  is  utilized  in  an  impedance  equation. 
The  equation  generates  an  impedance  matrix  so  that  the 
"method  of  moments"  can  be  used  to  solve  for  the  current 
distribution  or  the  input  impedance  of  an  antenna  mounted 
on  a  dielectric  substrate. 


II.   NATURE  OF  THE  PROBLEM 

Since  the  electrostatic  problem  has  been  solved,  a 
method  for  solving  the  dynamic  problem  is  sought  using 
image  theory.  The  dynamic  situation  is  chosen  to  be  a 
half-wave  dipole  antenna  on  a  dielectric  substrate  of 
arbitrary  thickness  and  relative  permittivity  (e  ).  A 
wire  constitutes  the  antenna  whose  radius  and  length  are 
variable. 

A.   ON  THE  IMAGE  COEFFICIENT  METHOD 

Consider  one  small  segment  of  the  wire  as  an  element 
of  charge.   Assume  the  charge  to  be  a  short  distance  from 
the  dielectric.   The  charge  element  will  have  many  lines 
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of  flux  in  every  radial  direction.   However,  in  considering 
a  single  line  of  flux,  some  of  it  penetrates  the  dielectric, 
while  some  of  it  is  reflected  at  the  interface.   The  amount 
reflected  is  proportional  to  the  reflection  coefficient  (K) . 
Similarly,  the  part  penetrating  the  dielectric  is  propor- 
tional to  (1-K). 

Naturally,  many  reflections  and  transmissions  occur 
due  to  this  single  flux  line.   In  turn,  it  produces  many 
image  points.   The  diagram  serves  to  illustrate: 
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(Figure  1. ) 


Each  of  the  three  regions,  I, II, III,  may.  have  different 
dielectric  constants.   However,  in  this  case,  regions  I  and 
III  have  e  =1  (free  space).   Region  II  has  an  arbitrary  e 
In  solving  for  the  Green's  Function,  each  region  has  a 
different  image  representation. 

(a).   Image  representation  for  left  region  -  I. 
(b).   Image  representation  for  center  region  -  II. 
(c).   Image  representation  for  right  region  -  III. 
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(Figure  2. ) 
B.   SOLVING  FOR  THE  GREEN'S  FUNCTION 

Specifically,  in  the  problem  at  hand,  the  wire  actually 
lies  on  the  dielectric  at  the  interface  of  regions  II  and 
III.   Region  II  is  the  desired  region  with  which  to  work 
because  the  accumulation  of  charge  is  assumed  to  be  greater 
on  the  dielectric  side  of  the  wire  than  on  the  free  space 
side.   The  image  representation  for  this  region  appears  as 


-PO-P)^   -p(l-P)q    (|-f)(j   PVp)<j  pVrti 

where  the  reflection  coefficient  is  now  p  c  1  ~€*-    in  order 
to  avoid  confusion  with  the  propagation  constant  (k). 
Adding  terms  to  simplify: 

(a)*     For  t    :      O-P) 


(b).  For't,-.   p2(i-p)-p(i-p)   --  (\-e)(?x-n 


(c).  For  %: 


(d).    For  ^3  ■• 
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-  In  terms  of  a  summation,  the  Green's  Function  for  the 
dielectric  region  is  expressed  ,  ,  .  . 
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Yet,  with  a  single  wire,  variation  occurs  in  only  the 
X-direction.  Y 
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where 


The  dielectric  Green's  Function  then  becomes 


t]t.  «exp(-jk/u -*')'  +  (2*D)Z) 

4*/ (*-*•')* 


C.   PHASE  CONSIDERATION 

In  formulating  the  Green's  Function  for  the  static 
problem,  the  phase  of  each  image  point  need  not  be  consid-' 
ered.   However,  in  the  dynamic  problem,  image  point  phase 
must  be  examined.   With  the  charge  element  a  distance  "a" 
from  the  dielectric  in  Figure  1,  each  image  point  has  a 
phase  term  associated  with  it.   This  means  that  H^d  could 
not  be  expressed  as  a  simple  summation. 

On  the  other  hand,  phase  terms  are  eliminated  if  the 
charge  element  is  adjacent  to  the  dielectric.   Since  the 
wire  is  actually  adjoined  to  the  dielectric,  every  image 
point  appears  as  though  it  were  immersed  in  the  dielectric 
region.   For  instance,  consider  image  point  "q, "  which 
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results  from  image  path  CI.   This  CI  comes  from  flux  line  A, 
which  is  reflected  into  B  and  then  C.   The  image  path  length 
(CI)  is  equal  to  the  addition  of  dielectric  path  lengths  A 
and  B  through  simple  geometry.   In  other  words ,  the ' image 
path  "CI"  appears  to  be  in  the  dielectric  medium  due  to  A 
and  B.   Then  image  point  "q-,"  also  appears  to  originate  in 
the  dielectric  for  the  same  reason.   This  is  true  for  every 
image  point.   Hence,  all  image  points  appear  to  be  in  the 
same  medium  with  equal  phase.   The  Green's  Function,  then, 
retains  its  simple  form. 

D.   THE  USE  OP  l|r. 

The  Green's  Function  (^d)  of  this  section  was  derived 
in  a  manner  similar  to  Silvester's  solution  [1]  using  the 
method  of  partial  images.  Silvester,  of  course,  solved  a 
static  problem.  Yet,  it  seems  logical  to  assume  that  his 
method  should  be  applicable  to  a  dynamic  problem  through 
the  reasons  stated. 

The  method  of  .moments  is  used  to  find  a  numerical 
solution  to  the  input  impedance  of  the  antenna.   This  method 
utilizes  the  Green's  Function  in  solving  for  the  impedance 
matrix. 


12 


III.   THE  IMPEDANCE  MATRIX 

In  solving  for  the  current  distribution  or  the  feed- 
point  impedance  of  an  antenna,  the  matrix  approach  in  the 
method  of  moments  is  a  good  approximation.   The  key  is  to 
properly"  load  the  impedance  matrix.   Once  this  is  accom- 
plished, assumed  values  of  voltage  fill  the  voltage  matrix 
and  the  problem  is  solved. 


[I]  =  [Z]_1  [V] 


The  impedance  matrix  is  obtained  through  the  impedance 
equation.   The  equation  is  derived  by  Harrington  [2]  and 
the  details  are  given  in  Appendix  A.   The  Green's  Function 
(^j)  from  the  previous  section  is  employed  in  the  equation, 


^(S,m)  "  t(n,m)  +  t^m)] 


A.   CALCULATION  OF  THE  GREEN'S  FUNCTION 

In  order  to  express  the  Green's  Function  in  matrix  form. 
M/-,  must  be  integrated  along  the  wire.   Holding  the  source 
point  (n)  fixed  while  the  observation  point  (m)  ranges  over 
every  segment  generates  the  n —  column  of  the  matrix.   Then 
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st 
moving  the  source  point  to  the  (n+1) —  segment  and,  again 

allowing  the  observation  point  (m)  to  roam  along  the  wire 

generates  the  next  column,  etc.   Finally,  the  (m  x  n) 

matrix  is  filled,  where  "m"  is  chosen  numerically  equal 

to  "n". 

In   terms   of  equations: 
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There  is  a  different  Green's  Function  for  each  case. 
The  derivations  are  detailed  in  Appendix  B  while  the  results 
are  shown  here. 

For  case  1.  (m=n) 
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For  case  II.  (m^n) 
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B.   CALCULATION  OF  THE  IMPEDANCE  MATRIX 

The  impedance  values  are  now  found  using  the  calculated 
Green's  Functions  in  the  impedance  equation. 

Z(m,i03   JUuAfin^Cn  m)    +  J_P)'(nim) 
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Where  Axn  "  Axn*  Ai-rn   since  the  length  of  each  segment 
is  chosen  to  be  equal. 

Three  different  situations  arise  in  calculating  elements 
of  the  impedance  matrix.   The  first  is  the  "self-impedance" 
term.   Whenever  the  observation  point  and  source  point 
coincide,  a  self-impedance  term  is  generated.   These  terms 
occupy  the  diagonal  of  the  matrix  and  are  all  numerically 
equal. 

The  second  situation  is  the  "mutual-adjacent"  term. 
These  occur  when  the  observation  point  is  adjacent  to  the 
source  point.   Due  to  the  fact  that  the  distance  between 
"m"  and  "n"  is  always  the  length  of  one  segment  (Ait\)> 
every  mutual-adjacent  term  is  equal. 

The  third  situation  happens  when  the  observation  and 
source  points  are  one  or  more  segments  apart  and  are  called 
"mutual  impedance"  terms.   They  complete  the  matrix  and 
vary  in  magnitude. 
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1 .      Self -Impedance 
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2.      Mutual-Adjacent    Impedance 
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3.      Mutual   Impedance 
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C.   PROGRAMMING 

The  impedance  equations  of  the  previous  section  are 
programmed  for  computer  usage.   Of  course,  similarities 
exist  throughout  the  equations  making  the  programming 
easier. 

Once  the  proper  parameters  and  constants  are  selected, 
the  impedance  matrix  is  filled  and  then  inverted.   The 
values  of  voltage  in  the  voltage  matrix  are  chosen  to  be 
zero  everywhere  except  at  the  feedpoint,  where  the  source 
is  one  volt. 


-1 
[Z]   [V]  =  [I] 


The  resulting  matrix  [I]  will  be  the  current  distribution 
on  the  wire.   The  middle  value  of  [I],  the  feedpoint  value, 
is  the  feedpoint  admittance  since  the  source  was  one  volt. 
The  computer  program  is  written  to  solve  for  various  values 
of  feedpoint,  or  input,  admittances. 
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IV.   DATA  EVALUATION 

After  programming  the  equations,  meaningful  data  is 
sought  for  analysis.   Several  parameters  are  available: 
the  radius,  the  length,  the  dielectric  constant,  the 
thickness,  and  the  frequency.   In  order  to  evaluate  the 
validity  of  the  equations,  a  comparison  is  made  against 
Harrington's  solutions.    Hence,  the  dielectric  constant 
in  Region  II  is  chosen  to  be  e  =1.   The  length  is  equal  to 
half  of  a  free  space  wavelength  at  f  =  3  GHz.   The  length- 
to-diameter  ratio,  L/2a,  is  set  at  7^.2,  while  the  frequency 
is  varied. 

The  results  are  plotted  and  compare  favorably  to 
Harrington's  graphs.   Figures  3  and  4  show  conductance  and 
susceptance,  respectively,  versus  the  length-to-wavelength 
ratios.   It  was  found  that  at  least  a  forty-one  by  forty-one 
matrix  was  necessary  to  achieve  accuracy. 

A.   EXPERIMENTAL  PROCEDURE 

Experimentally,  a  thin  copper  strip  was  attached  to  a 
one-eighth  inch  thick  dielectric  substrate  with  e  =16.   The 
dielectric  was  mounted  perpendicular  onto  the  ground  plane 
so  that  the  thin  strip  would  be  fed  as  a  quarter-wave 
monopole  vice  a  half-wave  dipole.   The  antenna  was  driven 


Harrington,  R.P.,  Field  Computation  by  Moment  Methods 
p.  72,  The  MacMillan  Company,  1968^ 
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through  50  ohm  miniature  coax  and  input  impedance  was 
measured  with  a  network  analyzer. 

The  experimental  values  of  input  conductance  and  sus- 
ceptance  are  plotted  along  with  values  calculated  using  the 
.impedance  equations.   Obviously,  Figures  5  and  6  illustrate 
a  significant  difference  between  experimental  and  calculated 
results.   There  is  an  amplitude  variation  and  a  difference 
in  resonant  frequencies. 
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V.   CONCLUSIONS 

Prom  the  free-space  curves  which  match  Harrington's 
curves,  it  seems  that  the  Green's  Function  and  the 
impedance  equations  are  valid.   Of  course,  in  this  instance, 
e  =1  so  that  the  reflection  coefficient  (p)  is  zero  and 
every  term  containing  (p)  is  eliminated.   The  part  remaining 
is  simply  a  free-space  Green's  Pucntion.   Hence,  the  solution 
is  valid  for  a  half-wave  dipole  antenna  in  free-space. 
Unfortunately,  the  same  is  not  true  when  er>l« 

A.  EXPERIMENTAL  DIFFERENCES 

Aside  from  the  normal  attenuation  losses  of  experimental 
measurements,  it  seems  the  only  other  inaccuracy  would  come 
from  the  adhesive  used  to  attach  the  copper  strip  to  the 
dielectric.   This  would  influence  the  effective  dielectric 
constant  (e   )  and,  perhaps,  the  amplitude  of  the  admit- 
tance.  Therefore,  the  emphasis  is  switched  to  the 
calculation  procedures. 

B.  PROBLEM  IN  THE  GREEN'S  FUNCTION 

The  manner  in  which  the  Green's  Function  was  chosen 
allows  some  error  to  exist.  The  accumulation  of  charge 
was  assumed  to  be  greater  on  the  dielectric  side  of  the 
wire  rather  than  the  free-space  side.  Because  of  this, 
Region  II  was  selected  in  accordance  with  the  method  of 
partial  images  to  solve  for  the  necessary  Green's  Function. 
In  the  process,  the  propagation  constant  is  dependent  on 
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the  relative  permittivity  of  the  dielectric. 


k=  I* 


2   it 


/U    A/IT 

However,  from  the  physical  standpoint,  as  the  limit  of 
the  dielectric  thickness  approaches  zero,  the  Green's 
Function  should  reduce  to  a  free-space  solution.   Yet,  y, 
from  section  II-B  will  never  reduce  to  a  free-space  function 
due  to  the  dependence  of  the  propagation  constant  on  e   . 
This  suggests  the  use  of  both  a  free-space  propagation 
constant  and  a  dielectric  propagation  constant.   But  there 
does  not  seem  to  be  any  apparent  method  in  combining  the  two 
propagation  constants  with  a  Green's  Function  for  Region  II. 

Further  investigation  led  to  a  Green's  Function  for 
Region  III  which  incorporated  both  propagation  constants. 
It  reduces  to  the  free-space  function  as  zero  dielectric 
thickness  is  approached.   However,  the  solution  leads  to 
negative  conductance. 

Finally,  since  the  method  of  partial  images  can  produce 
valid  solutions  for  static  and  free-space  dynamic  problems, 
it  would  seem  that  with  some  additional  changes  this  same 
method  should  solve  all  dynamic  problems.   Unfortunately, 
the  final  answer  lies  beyond  the  research  efforts  of  this 
paper.   The  method  of  partial  images,  however,  warrants  the 
need  for  further  investigation. 
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APPENDIX  A 
DERIVATION  OF  THE  IMPEDANCE  EQUATION 

The  impedance  equation  is  developed  in  a  manner  similar 
to  Harrington's  derivation.   Formulating  the  problem: 


Es=  -j<oA  -  v$ 


(A-0 


K 


-      -ikR 
Js_e aS 

4tr  R 


(a -a) 


€ 


(rs  _e aS 

4-n-R 


(a -3) 


(A -4) 


r\  x  E.      =  -n  x  E.*      on  S 


(A-S) 


2      *     a 


n     °     "        (Ftaurt    A-l) 
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The  wire  is  broken  into  many  segments  subject  to  the 
following  conditions: 

(1)  Current  flow  is  assumed  to  be  along  the  wire  axis. 

(2)  To  the  axial  component  of  E  at  the  surface,  apply 
boundary  condition  (A-5). 

With  these  assumptions,  equations  (A-l)  through  (A-4) 

become 


Ejj    z   -  j  cj  A  j>    -  j*J?_        on  S  (A-0 

hi 

7  r  -;JkR 

A   =  (i  \KH.e       afi  (a-t) 

L      4** 

$    -  j_V   <r(tt_e a£  (A-8) 

£   J  4frR 

AtflS 

cr  -  -j_  al_  (a-<?) 


where  "x  "  is  the  "length  variable"  along  the  axis  of  the 
wire. 
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Figure  A-l  illustrates  the  division  of  the  wire  into 

"n"  segments.   According  to  the  notation,  any  n —  segment 

_  + 

starts  at  n,  is  centered  at  n3  and  terminates  at  n.   The 

distance  between  any  n  and  n  is  denoted  by  Axn  5  the 

length  of  the  segment.   Summing  integrals  over  many  small 

segments  of  ACn  approximates  an  integral  over  the  entire 

length.   Equations  (A-6)  through  (A-9)  are  then  further 

approximated 


tjjM    *    iU>Aj,(nO    -     I  (m)  -    §  (m)  (A-10) 

AJJm 


A(m)  ~  \x\~*   'LmAK  \jL_-  a  A     (A-u) 


n 


j  Ah     UttR 


(m)    ~       1     \         0"(n)\    e 61  (A-\l) 


n  in 

(T(S)    ^    -J_       I  (tm)  -  X  (n)  (A- 13) 

^  AK 
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Similarly 

<rU)   «   -J_    X(n)  -  X  (n-Q  (A-14) 

iW  A  Jin 

substituting  equation    (A-13)    into    (A-12),    and   equation 
(A-l^)    into   an   equation   similar  to    (A-12)    for       S?(m), 


1  Cm)   =    -   }     V^    ICn+O-ICrQ    \     _C <dL  Q       (A -15) 


$(*)  =  -  i    \     Icn^-Icn-o   \    e        d&     (A-lO 


jO>t__j  A^a  J        +ffft 


it  is  often  convenient  to  express  the  integral  parts  of  the 
above  equations  in  the  form  of  Green's  Functions. 


t(nm)r  \        \    e"A    Ai 


n       )   .  4  «  R. 
Mr 


Where    Rm=  the  distance  from  the  midpoint  of  the  m — ,  , 
segment  to  the  integration  point  on  the  n — 


Using  this  definition,  equations  (A-ll),  (A-15)  and 
(A-16)  become 


Ac**)  =  ^Y"  I^M^M     ^A'17^ 


h 
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§ 


+ 

Cm) 


fc-iti 


icoe  /_,  L  J 

n 

Now  substituting  equations  (A-17),  (A-18)  and  (A-19)  into 

equation  (A-10)  and  noting  that  A|(m)=  A(m)  •   A  j m 


"Ejj  (tn)  -    -JGJp 


A  Jm 
XCn)Aln'A$m  ^(n.m) 


Ail 


m  Z 1 

n 


-J 

A? 


-  1 


m 


jui 


ICMO^CS.^-Icn)^^^) 


_    n 


*     I 


ju>£ 


!(«)%,«  -Itn-ONTCn.A) 


n 


J] 


(A-£o) 


Two  further  approximations  are  made 


and 


ICn-OYln,m)  -  I(n)tln>) 


(A-21) 
(A-**) 


With  these  two  substituted,  equation  (A-20)  is  written  as 


ELJmAJL  = 


dwP^ 


Iw  AL'ASL  ^taivO 


+  U    \        I(n) 


te^-^^-tcS.m)*^,^ 


.     (A- 23) 


-*   y 
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or 


Ej(r*)Af 


*A 


Z.(nfi,n)  I(|V) 


where 


n 


"-+(V*)  -!(S»  +  tu,*)l     (A-24) 

Equation  (A-2^)  is  the  impedance  equation  which  applies 
for  the  case  of  (m=n)  self  impedance,  as  well  as  for  (m^n) 
mutual  impedance. 

Since  "m"  roams  along  all  "n"  segments  of  the  wire, 
there  is  a  set  of  "n"  linear  equations.   In  matrix  form: 


Eto-Al. 


o 

o 


E(r>vAln 


o    o     o 


© 

o 


Zm  Zn2  °  °     °  Znn 


^_in 

X  (0 

Xcn 

o 

0 
o 

jLnn 

_X(n^_ 

or 

[V]  =  [Z]  [I] 

If  [V]  and  [Z]  are  known,  the  current  distribution  may  be 
found, 

[I]  =  [Z]~1[V] 
as  well  as  the  feedpoint  impedance. 
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APPENDIX  B 


DERIVATION  OF  THE  GREEN'S  FUNCTION 


There  is  a  distinct  Green's  Function  for  each  case: 
when  the  observation  point  coincides  with  the  source  point 
(m=n);  and  when  the  observation  point  and  source  do  not 
coincide  (m^n) . 

Case   I.    (m=n) 


t(n,m)  -    _l_    \(1-P)       €. 


-A](^-A-'t+  o? 


J 


d/3C' 
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P  g 


a.   The  first  term  must  be  integrated  over  one  segment 

T 


{P^^ry^ 


Using  a  Maclaurin  series  expansion  on  the  first  term; 
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d  rtl      -    (\  -  P) 
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33 


since  the  first  two  terms  of  this  Maclaurin  expansion  give 
reasonably  good  results, 


'-    0-P) 
4-irAlr, 


ZU  {d.m-/L')+]u»- £')*+<£ 


-  ijVdf 


assuming  A  5n  *\0a 


=    (I-P) 

4wy, 


ZJUhrV.jkAJn 


2  TT 


b.   Integration  is  now  performed  on  the  second  term  of 
the  Green's  Function. 
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Finally,  the  Green's  Function  for  Case  I,  (m=n)  is 


t(n,m)  '     (l-p) 
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Case  II.  (m^n) 


+  (n,m)  =  J. 


-  <  k  (/im-  An) 

0-p^  e    cU» 
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-  1 


(l-p) 
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For  the  crudest  approximation,  assume  that 
/  v^rYi-z^rO  +  (2iD)    remains  constant  over  the  integration, 
Then,  the  Green's  Function  for  (m^n)  becomes 

tu.m)  =  d-p^  e 
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C  APPENDIX  D.   COMPUTER  PROGRAM 

C 

C      CCMPUTER  APPLICATION  ON  THE  METHOD 

C      OF  IMAGES  USING  REGION  II 

C 

C 

C 

IMPLICIT  REAL*4  (K) 

IMPLICIT  COMPLEX*8  (C,Z) 

DIMENSION  Z(41,4l),  ZMP(41),  Y(41) 

REAL*4  CAMPS(41) 

REAL  V(41  1/2.0*0.0,  1.0,20*0.0/ 

COMPLEX  CABS 
C 
C      THE  THICKNESS  (D)  IN  METERS 

D=0. 003175 
C 
C      THE  DIELECTRIC  CONSTANT  (ER) 

ER=16.0 
C 
C      THE  RADIUS  (A)  IN  METERS 

A=0.C0O5O0 
C 

20  READ(5,30,END=250 )  DATA 
3C  F0RMAT(F5.3) 

NUM=41 

MID=(NUM+1 )/2 
C 
C      THE  LENGTH  OF  THE  WIRE  IN  METERS 

ELNGTH=0. 13493750 
C 

CLTA= EL NGTH/( FLOAT { NUM ) 1 

F=( (3E03)*DATA)/ELNGTH 

P  1  =  3. 14159 

W=2.0*PI*F 

U=PI*4.0*l.E-7 

E0=L.0/(36.0*PI  l*l.E-9 
C 
C 

C  EQUATION    CONSTANTS 

C 

RhC=(l .0-ER)/(l  .O  +  ER) 

B  =  C0 

K=W.*SQRT(U*EO*ER) 

A1=W*U*DLTA**2 

A2=( 1.0-RHO)/(2.0*PI) 

A3=(l .0/DLTA)*ALOG(DLTA/A) 

A4=1.0/(2.0*PI1*(1 ,0-RHQ)**2 

A5=2.0/ (W*EO*ER) 

A6=A2/2  .0 

A7=A4/2.0 
•     A8=A5/2.0 

A9=K/2.0 

A10=DLTA**2 

A11=4.0*D**2 

AK1=K*DLTA 

CD1=CEXP(CMPLX(B,-AKL) 1/DLTA 

AK2=2.0*AK1 

CC2=CEXP(CMPLX(B,-AK2 ) )/ (2 .0*DLTA ) 

C1  =  CMPLX(  B,A1) 

C2  =  A3-CMPLX(B,  A9) 
C 

C  SUMMATION    OF    THE     EXPONENTIAL    WITh    DELTA     (DLTA) 

C 

C3=CMPLX(B,B) 

R2=0.0 

DC    40    J=l, 10 

J  1  =  2.0* J- 1 

J2=J**2 

R1=RH0**J1 

R2=R2+R1 
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c 
c 
c 
c 


c 
c 
c 
c 


c 
c 
c 
c 


c 
c 
c 
c 


X1  =  SQRT  (A10+AU*J2) 
X2=K*X1 

C4=CEXP(CMPLX(B,-X2)) 
C5=(R1*C4)/X1 

C3=C3+C5 
4G    CONTINUE 


SUMMATION    OF    THE    EXPONENTIAL    WITH    2*DELTA     (2*DLTA) 

C6=CMPLX(B,B) 

R2  =  0.0 

DO    50    1=1, 10 

11=2.0*1-1 

12=1**2 

R4=RH0**I1 

R3=R3+R4 

X3=SQRT(4.0*A10+A11*I2 ) 

V Ar K  ^  X  ^ 

C8=CEXP(CMPLX(B,-X4)  ) 
C<5  =  (R4*C8) /X3 
C6=C6+C9 
50    CONTINUE 

SUMMATION    OF    THE    EXPONENTIAL    WITH    2*1*0 

C20=CMPLX(B,B) 
DC    60     1=1, 10 
12=2.0*1-1 
14=1**2 
P10=RH0**I3 
X2G=2.0*I*D 
X21  =  K*X20 

C21=CEXF(CMPLX(B,-X21 ) ) 
C22  =  (R10*C21)  /X20 
C20=C20+C22 
60    CONTINUE 


CALCULATION    OF    THE     IMPEDANCE    MATRIX 


DO    110 
DO    100 
L=IABS( 
IF(L.EQ 
IFtL.EQ 


M=1,NUM 

N  =  l  ,NUM 

M-N) 

.0)    GO    TO 

.1  )    GO    TO 


90 
80 


FOR    MUTUAL     IMPEDANCE 

LP=L+1 

LM=L-1 

AK3=K*L*DLTA 

CD3=CEXP(CMPLX(B,-AK3) )/(L*DLTA) 

AK4=K*LP*DLTA 

CD4=CEXP(CMPLX(B,-AK4) )/( LP*DLTA) 

AK5=K*LM*DLTA 

CD5=CEXP(CMPLX(B,-AK5) )/(LM*DLTA) 

C10=CMPLX(B,B) 

C13=CMPLX(B, B) 

C16=CMPLX(B,B) 

R5=0.0 

DC    70    J=l,10 

J11=2.0*J-1 

J22=J**2 

R6=RH0**J11 

R5=R5+R6 

A12=(L*DLTA)**2 

X5=SQRT(A12+A11*J22) 

X6=K*X5 

Cil=CEXP(CMPLX(B,-X6) ) 
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c 
c 
c 
c 


c 
c 
c 
c 


c 
c 

c 

c 
c 
c 
c 


C12=(R6*C11) /X5 
C10=C  10+C12 
A13=(LP*DLTA)**2 
X7  =  SQRT( A13+A11*J22  ) 
X  F  —  K  ^  X  7 

C14=CEXP(CMPLX(B,-X8) ) 
Ci5=(R6*C14) /X7 
C13=C13+C15 
A14=(LM*DLTA)**2 
X9=SQRT(  A14+A11*J22) 
X  1 0  =  K.  -^  X  9 

C17=CEXP.(CMPLX(  B,-X10  )  ) 
C18=(R6*C17)/X9 
C  16  =  C16+C18 
70    CONTINUE 

Z21=C1*(A6*CD3-A7*C10) 
Z22=(CMPLX(B,-A5 ) ) *( A6*CD3-A7*C  10) 
Z23=(CMPLX(B, A3)  )* ( A6*CD4-A7 *C1 3  ) 
Z24=(CMPLX(B,A8) )* ( A6*CD5-A7*C1 6 ) 
ZiM,N)=     Z21    +    Z22    +    Z23    +    Z24 
GC    TO    100 


FOR    MUTUAL-ADJACENT    IMPEDANCE 

80    Z'31=-C1*(A6*CD1-A7'*C3) 

Z32=(CMPLX(B,-A5)  ) * ( A6 *CD1 -A7* C3  ) 
Z33=(CMPLX(B,A8)  ) * ( A6*CD2-A7*C6  ) 
Z34=(CMPLX(B,A8 )  )* ( A2*C2-A7*C20  ) 
Z(M,N)     =    Z3L  +  Z32+Z33+Z34 
GO    TO     100 


FOR    SELF    IMPEDANCE 

90  Z.11=C1*(A2*C2"A7.*C20) 

Z12=(CMPLX(B,-A5) ) * ( A2*C2-A7*C20 ) 

Z13=(CMPLX(B,A5) )*(A6*CD1-A7*C3) 

Z(MiN)  =  Z11+Z12+Z13 
100  CONTINUE 
110  CONTINUE 

THE  MATRIX  INVERSION 

CALL  CM1N1(NUM,Z,NUM,DETERM) 

MATRIX    MULTIPLICATION    OF    Z-INVERSE    AND    V 

DC     130     1=1, NUM 

ZMPd  )  -(0.0,0.0) 

DC     120    J=1,NUM 

ZMP(  I)  =  ZMP(I  J+Z(  I, J)*V( J ) 

CONTINUE 

CONTINUE 

WPITE(6,140) 

FORMAT ( ' 1* ) 

WRITE(6,  150) 

FCPMAT( 15(/) 

WRITE(6il60)     D,ER,A 

FORMAT (  10X, 'THE  THICKNESS  D 
l'DIELECTRIC  CONSTANT  ER«,4X 
2'    A'  ,  18X,F  12.8,//) 

WRITE(6,170)     W 

FCRMAT(//,10X,' THE    RADIAN    FREQUENCY    W «  , 10X , E  1  3  .7 , / /  ) 

WRITE(6, 180) 

FORMAT! 10X, 'THE  COMPLEX 
3' AMPS' , /) 

WRITE(6,190)  ZMP(MID) 
90  FORMAT (42X, El  1.5, 3X, Ell 

DO  200  1=1, NUM 

C£MPS( I )  =  CABS  (ZMP(  I  )  ) 


120 
130 

140 

150 

160 


170 
130 


) 


,  15X 
,  F12  . 


,F12.8, //,  10X,' THE  ', 
3,//,  10 X,  'THE  RADIUS'. 


CURRENT  AT  THE  FEEDPOINT  IN  ' 


5,3X) 
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c 
c 
c 


200 
2L0 
220 
230 
240 

250 
26C 


CONTINUE 

WRITE(6,2L0) 

FORMAT (//, LOX, 'THE  MAGNITUDE  OF  THE  CURRENT',/) 

WRITE(6.220)  CAMPS(MID) 

F0RMAT(42X,ELL.  5) 

WRITE(6,230)     DATA 

FORMAK //,LOX,'FOR  THE  LENGTH/LAMDA  RATIO  «,5X,F6.4,/) 

WRITE(6,240)  F 


F0RMAT(//,10X,'THE  FREQUENCY 


L6X, 1PE1C.4, //) 


CARD  PRINT  OUT  -  DATA 
kRITE(7,25C)  ZMP(MID) 
FORMAT( ELL. 5) 

GO  TO  20 
CONTINUE 
STOP 
END 


FOR  DRAW  SUBROUTINE 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 


c 
c 
c 


c 
c 
c 


LO 
20 


SUBROUTINE    CMINL     (  N  ,  A ,  NDI  M,  DETE  RM  ) 


CMINL  IS  A  SUBROUTINE  WHICH  WILL  ACCEPT  A  SINGLE 
PRECISION  COMPLEX  MATRIX  AND  RETURNS  THE  INVERSE 
OF    THE     MATRIX    IN    ITS    PLACE. 


6C 
80 
85 


LOO 
LOS 


N 
A 


NDIM      - 


THE     ORDER    OF    THE    MATRIX    TO    BE     INVERTED 
COMPLEX    SINGLE    PRECISION    INPUT    MATRIX 
(DESTROYED).     THE     INVERSE    OF    A     IS 
RETURNED     IN    ITS    PLACE. 
THE    SIZE     TO    WHICH    A    IS    DIMENSIONED 


COMPLEX    A (NDI M,NDIf)  , PIVOT ( LOO)  , AMAX , T , SWA P  ,DETERM , U 
INTEGER**    IPIVOK  LOO) ,  INDEX ( LOO, 2  ) 
REAL    TEMP,ALPHA( LOO) 

INITIALIZATION 

DETERM=CMPLX(  L. 0,0.0) 

DC    20    J  =  L,N 

ALPHA     (J)=0.ODO 

DO    LO    I  =  L,N 

ALPHA (J )  =  ALP  HA (J )+A(J,  I 

ALPHA( J)=SCRT(ALPHA(J) ) 

IPIVOK  J)  =  0 

DC    600     1  =  1,  N 


)*CONJG( A(J,I  )) 


SEARCH    FOR    PIVOT    ELEMENT 

AMAX=CMPLX(0. 0,0.0) 

DO     105    J=l,N 

IF    (  IPIVOKJ  )-L  )     60,L05,60 

DC    LOO    K=L,N 

IF    ( IPI VOT(K )-L)     80, 

TEMP=AMAX*CONJG(AMAX  ) 

IF(TEMP)     85,    85,     LOO 

IROW=J 

ICCLUM=K 

AMAX=A( J,K) 

CONTINUE 

CONTINUE 

IPIVOK  ICOLUM) 


100,740 

A(J,K)*CONJG(A(J,K) ) 


IPIVOK ICOLUM)+L 


INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

IF(  IROW-ICOLUM)  L40,  260,  L40 
L4C  DETERM=DETERM 
DO  200  L=L,N 
SWAP=A( IROW,L) 


40 


c 
c 
c 


c 
c 
c 


c 
c 
c 


AC  IROW,L)=A(  ICOLUM, L) 
200  AC  IC0LUM,L)  =  SWAP 

SkAP=ALPHA(  IROW) 

ALPHA ( IROW)=ALPHA( ICOLUM) 

ALPHA( ICOLUM)=SWAP 
260  INDEXC I  ,1) =IROW 

INDEX( I,2)=IC0LUM 

PIVOT  CI  )  =  A( ICOLUM,  ICOLUM) 

U=PIVOT(I) 

TEMP=PIVOT(I ) *CONJG(PIVOT( 

IF(TEMP)  330,  720,  330 


I  )  ) 


DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 

330  AC ICOLUM, ICOLUM) =CNPLX  (1  .0,0.0  ) 

DO  350  L=1,N 

U=PIVOT  (I  ) 
350  AC ICOLUM, L ) = A ( I COLUM , L ) /U 

REDUCE  NON-PIVOT  ROWS 

380  DO  550  L1=1,N 

IFCL1-IC0LUM)  400,  550,  400 
400  T=A(L1, ICOLUM) 

ACL1, ICOLUM )=CMPLX(C.  0,0.0) 

DC  450  L=lfN 

U=A( ICOLUM, L) 
450  A(LI,L)=A(L1,L)-U*T 
550  CCNTINUE 
60C  CONTINUE 

INTERCHANGE  COLUMNS 

620  DO  710  1=1, N 
L=N+1-I 

IFC INDEXCL ,L)-INDEX(L,2) )  630,  7L0,  630 
63C  JRCW=INDEX(L,  I) 

JCGLUM=INDEX(L,2 ) 

DO  705  K=1,N 

SWAP=A(K, JROW) 

A(K,JR0W) =A(K, JCOLUM) 

A(K,JCOLUM)=SWAP 
705  CCNTINUE 
7  1C  CCNTINUE 

RETURN 
720  WPITE(6,730) 

73C  F0RMATC20H   MATRIX  IS  SINGULAR) 
740  RETURN 

END 
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